from cmath import isinf
from distutils import dist
import math
import networkx as nx
from numba import cuda
import numpy as np
from . import floyd_in_parallel


def calculate_fstate_shortest_path_without_gs_relaying(
        output_dynamic_state_dir,
        time_since_epoch_ns,
        num_satellites,
        num_ground_stations,
        sat_net_graph_only_satellites_with_isls,
        num_isls_per_sat,
        gid_to_sat_gsl_if_idx,
        ground_station_satellites_in_range_candidates,
        sat_neighbor_to_if,
        prev_fstate,
        enable_verbose_logs
):

    if cuda.is_available():
        free=cuda.current_context().get_memory_info()[0]
        if free>=1000000000:
            return calculate_fstate_shortest_path_without_gs_relaying_gpu(
                output_dynamic_state_dir,
                time_since_epoch_ns,
                num_satellites,
                num_ground_stations,
                sat_net_graph_only_satellites_with_isls,
                num_isls_per_sat,
                gid_to_sat_gsl_if_idx,
                ground_station_satellites_in_range_candidates,
                sat_neighbor_to_if,
                prev_fstate,
                enable_verbose_logs
            )

    # Calculate shortest path distances
    if enable_verbose_logs:
        print("  > Calculating Floyd-Warshall for graph without ground-station relays")
    # (Note: Numpy has a deprecation warning here because of how networkx uses matrices)

    dist_sat_net_without_gs = nx.floyd_warshall_numpy(sat_net_graph_only_satellites_with_isls)

    # Forwarding state
    fstate = {}

    # Now write state to file for complete graph
    output_filename = output_dynamic_state_dir + "/fstate_" + str(time_since_epoch_ns) + ".txt"
    if enable_verbose_logs:
        print("  > Writing forwarding state to: " + output_filename)
    with open(output_filename, "w+") as f_out:

        # Satellites to ground stations and ground stations to satellites
        # From the satellites attached to the destination ground station,
        # select the one which promises the shortest path to the destination ground station (getting there + last hop)
        dist_satellite_to_ground_station = {}
        for curr in range(num_satellites):
            for dst_gid in range(num_ground_stations):
                dst_gs_node_id = num_satellites + dst_gid

                # Among the satellites in range of the destination ground station,
                # find the one which promises the shortest distance
                possible_dst_sats = ground_station_satellites_in_range_candidates[dst_gid]
                possibilities = []
                for b in possible_dst_sats:
                    if not math.isinf(dist_sat_net_without_gs[(curr, b[1])]):  # Must be reachable
                        possibilities.append(
                            (
                                dist_sat_net_without_gs[(curr, b[1])] + b[0],
                                b[1]
                            )
                        )
                possibilities = list(sorted(possibilities))

                # By default, if there is no satellite in range for the
                # destination ground station, it will be dropped (indicated by -1)
                sat_to_gs_next_hop_decision = (-1, -1, -1)
                gs_to_sat_next_hop_decision = (-1, -1, -1)
                distance_to_ground_station_m = float("inf")
                if len(possibilities) > 0:
                    dst_sat = possibilities[0][1]
                    distance_to_ground_station_m = possibilities[0][0]
                    gs_to_sat_next_hop_decision = (
                        dst_sat,
                        0,
                        num_isls_per_sat[curr]+gid_to_sat_gsl_if_idx[dst_gid]
                    )

                    # If the current node is not that satellite, determine how to get to the satellite
                    if curr != dst_sat:

                        # Among its neighbors, find the one which promises the
                        # lowest distance to reach the destination satellite
                        best_distance_m = 1000000000000000
                        for neighbor_id in sat_net_graph_only_satellites_with_isls.neighbors(curr):
                            distance_m = (
                                    sat_net_graph_only_satellites_with_isls.edges[(curr, neighbor_id)]["weight"]
                                    +
                                    dist_sat_net_without_gs[(neighbor_id, dst_sat)]
                            )
                            if distance_m < best_distance_m:
                                sat_to_gs_next_hop_decision = (
                                    neighbor_id,
                                    sat_neighbor_to_if[(curr, neighbor_id)],
                                    sat_neighbor_to_if[(neighbor_id, curr)]
                                )
                                best_distance_m = distance_m

                    else:
                        # This is the destination satellite, as such the next hop is the ground station itself
                        sat_to_gs_next_hop_decision = (
                            dst_gs_node_id,
                            num_isls_per_sat[dst_sat] + gid_to_sat_gsl_if_idx[dst_gid],
                            0
                        )

                # In any case, save the distance of the satellite to the ground station to re-use
                # when we calculate ground station to ground station forwarding
                dist_satellite_to_ground_station[(curr, dst_gs_node_id)] = distance_to_ground_station_m

                # Write to forwarding state
                if not prev_fstate or prev_fstate[(curr, dst_gs_node_id)] != sat_to_gs_next_hop_decision:
                    f_out.write("%d,%d,%d,%d,%d\n" % (
                        curr,
                        dst_gs_node_id,
                        sat_to_gs_next_hop_decision[0],
                        sat_to_gs_next_hop_decision[1],
                        sat_to_gs_next_hop_decision[2]
                    ))

                if not prev_fstate or prev_fstate[(dst_gs_node_id,curr)] != gs_to_sat_next_hop_decision:
                    f_out.write("%d,%d,%d,%d,%d\n" % (
                        dst_gs_node_id,
                        curr,
                        gs_to_sat_next_hop_decision[0],
                        gs_to_sat_next_hop_decision[1],
                        gs_to_sat_next_hop_decision[2]
                    ))
                
                fstate[(curr, dst_gs_node_id)] = sat_to_gs_next_hop_decision
                fstate[(dst_gs_node_id, curr)] = gs_to_sat_next_hop_decision

        # Ground stations to ground stations
        # Choose the source satellite which promises the shortest path
        for src_gid in range(num_ground_stations):
            for dst_gid in range(num_ground_stations):
                if src_gid != dst_gid:
                    src_gs_node_id = num_satellites + src_gid
                    dst_gs_node_id = num_satellites + dst_gid

                    # Among the satellites in range of the source ground station,
                    # find the one which promises the shortest distance
                    possible_src_sats = ground_station_satellites_in_range_candidates[src_gid]
                    possibilities = []
                    for a in possible_src_sats:
                        best_distance_offered_m = dist_satellite_to_ground_station[(a[1], dst_gs_node_id)]
                        if not math.isinf(best_distance_offered_m):
                            possibilities.append(
                                (
                                    a[0] + best_distance_offered_m,
                                    a[1]
                                )
                            )
                    possibilities = sorted(possibilities)

                    # By default, if there is no satellite in range for one of the
                    # ground stations, it will be dropped (indicated by -1)
                    next_hop_decision = (-1, -1, -1)
                    if len(possibilities) > 0:
                        src_sat_id = possibilities[0][1]
                        next_hop_decision = (
                            src_sat_id,
                            0,
                            num_isls_per_sat[src_sat_id] + gid_to_sat_gsl_if_idx[src_gid]
                        )

                    # Update forwarding state
                    if not prev_fstate or prev_fstate[(src_gs_node_id, dst_gs_node_id)] != next_hop_decision:
                        f_out.write("%d,%d,%d,%d,%d\n" % (
                            src_gs_node_id,
                            dst_gs_node_id,
                            next_hop_decision[0],
                            next_hop_decision[1],
                            next_hop_decision[2]
                        ))
                    fstate[(src_gs_node_id, dst_gs_node_id)] = next_hop_decision
        
        ## Satellites to satellites
        ## From satellites attach to the source satellites,
        ## choose the one which promises the shortest path
        for src_sat_id in range(num_satellites):
            for dst_sat_id in range(num_satellites):

                next_hop_decision=(-1,-1,-1)
                if src_sat_id!=dst_sat_id:
                    distance_min=math.inf
                    for src_neighbor_id in sat_net_graph_only_satellites_with_isls.neighbors(src_sat_id):
                        distance=(
                            dist_sat_net_without_gs[(src_neighbor_id,dst_sat_id)]
                          +
                            sat_net_graph_only_satellites_with_isls.edges[(src_sat_id,src_neighbor_id)]["weight"]
                        )
                        if distance<distance_min:
                            distance_min=distance
                            next_hop_decision=(
                                src_neighbor_id,
                                sat_neighbor_to_if[(src_sat_id,src_neighbor_id)],
                                sat_neighbor_to_if[(src_neighbor_id,src_sat_id)]
                            )
                if not prev_fstate or prev_fstate[(src_sat_id,dst_sat_id)] != next_hop_decision:
                    f_out.write("%d,%d,%d,%d,%d\n" %(
                        src_sat_id,
                        dst_sat_id,
                        next_hop_decision[0],
                        next_hop_decision[1],
                        next_hop_decision[2]
                    ))

                fstate[(src_sat_id,dst_sat_id)] = next_hop_decision


    # Finally return result
    return fstate


def calculate_fstate_shortest_path_with_gs_relaying(
        output_dynamic_state_dir,
        time_since_epoch_ns,
        num_satellites,
        num_ground_stations,
        sat_net_graph,
        num_isls_per_sat,
        gid_to_sat_gsl_if_idx,
        sat_neighbor_to_if,
        prev_fstate,
        enable_verbose_logs
):
    if cuda.is_available():
        free=cuda.current_context().get_memory_info()[0]
        if free>=1000000000:
            return calculate_fstate_shortest_path_with_gs_relaying_gpu(
                output_dynamic_state_dir,
                time_since_epoch_ns,
                num_satellites,
                num_ground_stations,
                sat_net_graph,
                num_isls_per_sat,
                gid_to_sat_gsl_if_idx,
                sat_neighbor_to_if,
                prev_fstate,
                enable_verbose_logs
            )

    # Calculate shortest paths
    if enable_verbose_logs:
        print("  > Calculating Floyd-Warshall for graph including ground-station relays")
    # (Note: Numpy has a deprecation warning here because of how networkx uses matrices)
    dist_sat_net = nx.floyd_warshall_numpy(sat_net_graph)

    # Forwarding state
    fstate = {}

    # Now write state to file for complete graph
    output_filename = output_dynamic_state_dir + "/fstate_" + str(time_since_epoch_ns) + ".txt"
    if enable_verbose_logs:
        print("  > Writing forwarding state to: " + output_filename)
    with open(output_filename, "w+") as f_out:

        # Satellites and ground stations to ground stations
        for current_node_id in range(num_satellites + num_ground_stations):
            for dst_node_id in range(num_satellites + num_ground_stations):

                # Cannot forward to itself
                if current_node_id != dst_node_id:

                    # Among its neighbors, find the one which promises the
                    # lowest distance to reach the destination satellite
                    next_hop_decision = (-1, -1, -1)
                    best_distance_m = 1000000000000000
                    for neighbor_id in sat_net_graph.neighbors(current_node_id):

                        # Any neighbor must be reachable
                        if math.isinf(dist_sat_net[(current_node_id, neighbor_id)]):
                            raise ValueError("Neighbor cannot be unreachable")

                        # Calculate distance = next-hop + distance the next hop node promises
                        distance_m = (
                            sat_net_graph.edges[(current_node_id, neighbor_id)]["weight"]
                            +
                            dist_sat_net[(neighbor_id, dst_node_id)]
                        )
                        if (
                                not math.isinf(dist_sat_net[(neighbor_id, dst_node_id)])
                                and
                                distance_m < best_distance_m
                        ):

                            # Check node identifiers to determine what are the
                            # correct interface identifiers
                            if current_node_id >= num_satellites and neighbor_id < num_satellites:  # GS to sat.
                                my_if = 0
                                next_hop_if = (
                                    num_isls_per_sat[neighbor_id]
                                    +
                                    gid_to_sat_gsl_if_idx[current_node_id - num_satellites]
                                )

                            elif current_node_id < num_satellites and neighbor_id >= num_satellites:  # Sat. to GS
                                my_if = (
                                    num_isls_per_sat[current_node_id]
                                    +
                                    gid_to_sat_gsl_if_idx[neighbor_id - num_satellites]
                                )
                                next_hop_if = 0

                            elif current_node_id < num_satellites and neighbor_id < num_satellites:  # Sat. to sat.
                                my_if = sat_neighbor_to_if[(current_node_id, neighbor_id)]
                                next_hop_if = sat_neighbor_to_if[(neighbor_id, current_node_id)]

                            else:  # GS to GS
                                raise ValueError("GS-to-GS link cannot exist")

                            # Write the next-hop decision
                            next_hop_decision = (
                                neighbor_id,  # Next-hop node identifier
                                my_if,        # My outgoing interface id
                                next_hop_if   # Next-hop incoming interface id
                            )

                            # Update best distance found
                            best_distance_m = distance_m

                    # Write to forwarding state
                    if not prev_fstate or prev_fstate[(current_node_id, dst_node_id)] != next_hop_decision:
                        f_out.write("%d,%d,%d,%d,%d\n" % (
                            current_node_id,
                            dst_node_id,
                            next_hop_decision[0],
                            next_hop_decision[1],
                            next_hop_decision[2]
                        ))
                    fstate[(current_node_id, dst_node_id)] = next_hop_decision

    # Finally return result
    return fstate

@cuda.jit
def calculate_satellites_to_satellites(
    num_satellites,
    adjacency_matrix,
    dist_sat_net,
    sat_neighbor_index_list,
    sat_neighbor_array_list,
    sat_to_sat_if,
    next_hop_decision_matrix
):
    src_sat_id,dst_sat_id=cuda.grid(2)
    if src_sat_id>=num_satellites or dst_sat_id>=num_satellites:
        return
    
    next_hop_decision=(-1,-1,-1)
    index_begin= 0 if src_sat_id==0 else sat_neighbor_index_list[src_sat_id-1]
    index_end=sat_neighbor_index_list[src_sat_id]

    distance_min=math.inf
    for i in range(index_begin,index_end):
        sat_neighbor=sat_neighbor_array_list[i]
        distance=adjacency_matrix[src_sat_id,sat_neighbor]+dist_sat_net[sat_neighbor,dst_sat_id]
        if distance<distance_min:
            distance_min=distance
            next_hop_decision=(
                sat_neighbor,
                sat_to_sat_if[src_sat_id,sat_neighbor],
                sat_to_sat_if[sat_neighbor,src_sat_id]
            )
    
    for i in range(3):
        next_hop_decision_matrix[src_sat_id,dst_sat_id][i]=next_hop_decision[i]

@cuda.jit
def calculate_path_between_satellites_and_ground_station(
    num_satellites,
    num_ground_stations,
    gs_neighbor_index_list,
    gs_neighbor_id_list,
    gs_neighbor_distance_list,
    dist_sat_net,
    gid_to_sat_gsl_idx,
    num_isls_per_sat,
    dist_sat_to_gs,
    next_hop_decision_matrix
):
    sat_id,gs_id=cuda.grid(2)
    if sat_id>=num_satellites or gs_id>=num_ground_stations:
        return
    
    gs_node_id=gs_id+num_satellites
    index_start= 0 if gs_id==0 else gs_neighbor_index_list[gs_id-1]
    index_end=gs_neighbor_index_list[gs_id]

    distance_min=math.inf
    dst_sat=-1
    for i in range(index_start,index_end):
        gs_neighbor=gs_neighbor_id_list[i]
        distance=gs_neighbor_distance_list[i]+dist_sat_net[sat_id,gs_neighbor]

        if distance<distance_min:
            distance_min=distance
            dst_sat=gs_neighbor
    
    ## calculate next hop
    sat_to_gs_next_hop_decision=(-1,-1,-1)
    gs_to_sat_next_hop_decision=(-1,-1,-1)
    if dst_sat!=-1:
        gs_to_sat_next_hop_decision=(
            gs_neighbor,
            0,
            num_isls_per_sat[gs_neighbor]+gid_to_sat_gsl_idx[gs_id]
        )
        if dst_sat!=sat_id:
            sat_to_gs_next_hop_decision=(
                next_hop_decision_matrix[sat_id,gs_neighbor][0],
                next_hop_decision_matrix[sat_id,gs_neighbor][1],
                next_hop_decision_matrix[sat_id,gs_neighbor][2]
            )
        else:
            sat_to_gs_next_hop_decision=(
                gs_node_id,
                num_isls_per_sat[dst_sat]+gid_to_sat_gsl_idx[gs_id],
                0
            )
    
    dist_sat_to_gs[sat_id,gs_id]=distance_min
    for i in range(3):
        next_hop_decision_matrix[sat_id,gs_node_id][i]=sat_to_gs_next_hop_decision[i]
        next_hop_decision_matrix[gs_node_id,sat_id][i]=gs_to_sat_next_hop_decision[i]

@cuda.jit
def calculate_ground_station_to_ground_station(
    num_satellites,
    num_ground_stations,
    gs_neighbor_index_list,
    gs_neighbor_id_list,
    gs_neighbor_distance_list,
    dist_sat_to_gs,
    next_hop_decision_matrix
):
    src_gs_id,dst_gs_id=cuda.grid(2)
    if src_gs_id>=num_ground_stations or dst_gs_id>=num_ground_stations:
        return 

    src_node_id=src_gs_id+num_satellites
    dst_node_id=dst_gs_id+num_satellites

    index_start=0 if dst_gs_id ==0 else gs_neighbor_index_list[dst_gs_id]
    index_end=gs_neighbor_index_list[dst_gs_id]

    distance_min=math.inf
    dst_sat=-1
    for i in range(index_start,index_end):
        gs_neighbor=gs_neighbor_id_list[i]
        distance=gs_neighbor_distance_list[i]+dist_sat_to_gs[gs_neighbor,src_gs_id]
        if distance_min>distance:
            distance_min=distance
            dst_sat=gs_neighbor
    
    next_hop_decision=(-1,-1,-1)
    if dst_sat!=-1:
        next_hop_decision=(
            next_hop_decision_matrix[src_node_id,dst_sat][0],
            next_hop_decision_matrix[src_node_id,dst_sat][1],
            next_hop_decision_matrix[src_node_id,dst_sat][2]
        )
    
    for i in range(3):
        next_hop_decision_matrix[src_node_id,dst_node_id][i]=next_hop_decision[i]

## calculate shortest path through CUDA
def calculate_fstate_shortest_path_without_gs_relaying_gpu(
    output_dynamic_state_dir,
    time_since_epoch_ns,
    num_satellites,
    num_ground_stations,
    sat_net_graph_only_satellites_with_isls,
    num_isls_per_sat,
    gid_to_sat_gsl_if_idx,
    ground_station_satellites_in_range_candidates,
    sat_neighbor_to_if,
    prev_fstate,
    enable_verbose_logs
):
     # Calculate shortest path distances
    if enable_verbose_logs:
        print("  > Calculating Floyd-Warshall for graph without ground-station relays")
    # (Note: Numpy has a deprecation warning here because of how networkx uses matrices)

    ## get adjacency matrix
    adjacency_matrix=nx.to_numpy_array(sat_net_graph_only_satellites_with_isls,nonedge=np.inf,dtype=np.float32)
    np.fill_diagonal(adjacency_matrix,0)
    dist_sat_net_without_gs_dev=floyd_in_parallel.floyd_gpu(adjacency_matrix)

    ## neighbor structure
    sat_neighbor_index_list=np.cumsum(sum(adjacency_matrix!=np.inf))
    sat_neighbor_array_list=np.array(np.where(adjacency_matrix!=np.inf))

    ## ground station neighbor
    gs_neighbor_index_list=np.cumsum(np.array([len(x) for x in ground_station_satellites_in_range_candidates]))

    ## convenge data
    sat_neighbor_length=sat_neighbor_index_list[-1]
    gs_neighbor_length=gs_neighbor_index_list[-1]
    
    sum_int=num_satellites+sat_neighbor_length+num_satellites*num_satellites+num_ground_stations+gs_neighbor_length+num_ground_stations+num_satellites*num_satellites+(num_satellites+num_ground_stations)*(num_satellites+num_ground_stations)*3
    sum_float=num_satellites*num_satellites+gs_neighbor_length+num_satellites*num_ground_stations

    data_int=np.full(sum_int,-1,dtype=np.int32)
    data_float=np.full(sum_float,-1,dtype=np.float64)
    
    index_start=0
    for i in range(0,num_satellites):
        data_int[index_start+i]=sat_neighbor_index_list[i]
    
    index_start=index_start+num_satellites
    for i in range(0,sat_neighbor_length):
        data_int[index_start+i]=sat_neighbor_array_list[1][i]
    
    index_start=index_start+sat_neighbor_length
    sat_to_sat_if=sat_neighbor_to_if.reshape(num_satellites*num_satellites)
    for i in range(0,num_satellites*num_satellites):
        data_int[index_start+i]=sat_to_sat_if[i]
    
    index_start=index_start+num_satellites*num_satellites
    for i in range(num_ground_stations):
        data_int[index_start+i]=gs_neighbor_index_list[i]
    
    index_start=index_start+gs_neighbor_length
    keys=range(0,num_ground_stations)
    for i in keys:
        for j in ground_station_satellites_in_range_candidates[i]:
            data_int[index_start]=j[1]
            index_start=index_start+1
    
    for i in range(num_ground_stations):
        data_int[index_start]=gid_to_sat_gsl_if_idx[i]

    for i in range(num_satellites):
        data_int[index_start+i]=num_isls_per_sat[i]
    
    index_start=index_start+num_satellites

    index_start=0
    for i in range(num_satellites):
        for j in range(num_satellites):
            data_float[index_start]=adjacency_matrix[i][j]
            index_start=0
    
    for i in range(num_ground_stations):
        for j in ground_station_satellites_in_range_candidates[i]:
            data_float[index_start]=j[0]
            index_start=index_start+1
    
    data_int_dev=cuda.to_device(data_int)
    data_float_dev=cuda.to_device(data_float)

    ## extract device array
    index_start=0
    sat_neighbor_index_dev=data_int_dev[index_start:index_start+num_satellites]
    index_start=index_start+num_satellites
    sat_neighbor_id_dev=data_int_dev[index_start:index_start+sat_neighbor_length]
    index_start=index_start+sat_neighbor_length
    sat_to_sat_if_dev=data_int_dev[index_start:index_start+num_satellites*num_satellites].reshape(num_satellites,num_satellites)
    index_start=index_start+num_satellites*num_satellites
    gs_neighbor_index_dev=data_int_dev[index_start:index_start+num_ground_stations]
    index_start=index_start+num_ground_stations
    gs_neighbor_id_dev=data_int_dev[index_start:index_start+gs_neighbor_length]
    index_start=index_start+gs_neighbor_length
    gid_to_sat_gsl_if_idx_dev=data_int_dev[index_start:index_start+num_ground_stations]
    index_start=index_start+num_ground_stations
    num_isls_per_sat_dev=data_int_dev[index_start:index_start+num_satellites]
    index_start=index_start+num_satellites
    next_hop_decision_dev=data_int_dev[index_start:index_start+((num_satellites+num_ground_stations)**2)*3].reshape((num_satellites+num_ground_stations),(num_satellites+num_ground_stations),3)

    index_start=0
    adjacency_matrix_dev=data_float_dev[index_start:index_start+num_satellites**2].reshape(num_satellites,num_satellites)
    index_start=index_start+num_satellites**2
    gs_neighbor_distance_dev=data_float_dev[index_start:index_start+gs_neighbor_length]
    index_start=index_start+gs_neighbor_length
    dist_sat_to_gs_dev=data_float_dev[index_start:index_start+num_satellites*num_ground_stations].reshape(num_satellites,num_ground_stations)

    ## invoke kernel
    TPB=(32,32)
    ## calculate the shortest path between satellites
    col=int((num_satellites+31)/32)
    row=int((num_satellites+31)/32)
    BPG=(row,col)
    cuda.synchronize()
    calculate_satellites_to_satellites[BPG,TPB](
        num_satellites,
        adjacency_matrix_dev,
        dist_sat_net_without_gs_dev,
        sat_neighbor_index_dev,
        sat_neighbor_id_dev,
        sat_to_sat_if_dev,
        next_hop_decision_dev
    )

    cuda.synchronize()
    row=int((num_satellites+31)/32)
    col=int((num_ground_stations+31)/32)
    BPG=(row,col)
    calculate_path_between_satellites_and_ground_station[BPG,TPB](
        num_satellites,
        num_ground_stations,
        gs_neighbor_index_dev,
        gs_neighbor_id_dev,
        gs_neighbor_distance_dev,
        dist_sat_net_without_gs_dev,
        gid_to_sat_gsl_if_idx_dev,
        num_isls_per_sat_dev,
        dist_sat_to_gs_dev,
        next_hop_decision_dev
    )

    cuda.synchronize()
    row=int((num_ground_stations+31)/32)
    col=int((num_ground_stations+31)/32)
    BPG=(row,col)
    calculate_ground_station_to_ground_station[BPG,TPB](
        num_satellites,
        num_ground_stations,
        gs_neighbor_index_dev,
        gs_neighbor_id_dev,
        gs_neighbor_distance_dev,
        dist_sat_to_gs_dev,
        next_hop_decision_dev
    )
    cuda.synchronize()

    next_hop_decision_matrix=next_hop_decision_dev.copy_to_host()
    fstate={}
    output_filename = output_dynamic_state_dir + "/fstate_" + str(time_since_epoch_ns) + ".txt"
    with open(output_filename,"w+") as f_out:
        for src in range(num_satellites+num_ground_stations):
            for dst in range(num_satellites+num_ground_stations):
                next_hop_decision=(
                    next_hop_decision_matrix[src,dst][0],
                    next_hop_decision_matrix[src,dst][1],
                    next_hop_decision_matrix[src,dst][2]
                )
                if not prev_fstate or prev_fstate[(src,dst)] != next_hop_decision:
                    f_out.write("%d,%d,%d,%d,%d\n" % (
                        src,
                        dst,
                        next_hop_decision[0],
                        next_hop_decision[1],
                        next_hop_decision[2],
                    ))
                fstate[(src,dst)]=next_hop_decision

    return fstate

@cuda.jit
def calculate_shortest_path_only_gsl_between_node(
    num_satellites,
    num_ground_stations,
    node_neighbor_index_list,
    node_neighbor_array_list,
    adjacency_matrix,
    dist_node_to_node,
    sat_to_sat_if,
    gid_to_sat_gsl_if_idx,
    num_isls_per_sat,
    next_hop_decision_matrix
):
    src_node_id,dst_node_id=cuda.grid(2)
    num_node=num_satellites+num_ground_stations
    if src_node_id>=num_node or dst_node_id>=num_node:
        return

    index_begin=0 if src_node_id==0 else node_neighbor_index_list[src_node_id-1]
    index_end=node_neighbor_index_list[src_node_id]

    next_hop_decision=(-1,-1,-1)
    if not math.isinf(adjacency_matrix[src_node_id,dst_node_id]):
        src_neighbor=-1
        distance_min=math.inf
        for i in range(index_begin,index_end):
            neighbor=node_neighbor_array_list[i]
            distance=adjacency_matrix[src_node_id,src_neighbor]+dist_node_to_node[src_neighbor,dst_node_id]
            if distance<distance_min:
                src_neighbor=neighbor
                distance_min=distance
    
        if src_neighbor != -1:
            ## sat to sat
            if src_node_id<num_satellites and src_neighbor>=num_satellites:
                curr_if=sat_to_sat_if[src_node_id,src_node_id]
                next_if=sat_to_sat_if[src_neighbor,src_node_id]
            
            ## sat to gs
            elif src_node_id<num_satellites and src_neighbor>=num_satellites:
                curr_if=num_isls_per_sat[src_node_id]+gid_to_sat_gsl_if_idx[src_neighbor-num_satellites]
                next_if=0

            ## gs to sat
            elif src_node_id>=num_satellites and src_neighbor<num_satellites:
                curr_if=0
                next_if=num_isls_per_sat[src_node_id]+gid_to_sat_gsl_if_idx[src_neighbor-num_satellites]
            
            ## gs to gs
            else:
                raise ValueError("GS-to-GS link cannot exist")

            next_hop_decision=(
                src_neighbor,
                curr_if,
                next_if
            )

            for i in range(3):
                next_hop_decision_matrix[src_node_id,dst_node_id][i]=next_hop_decision[i]

def calculate_fstate_shortest_path_with_gs_relaying_gpu(
        output_dynamic_state_dir,
        time_since_epoch_ns,
        num_satellites,
        num_ground_stations,
        sat_net_graph,
        num_isls_per_sat,
        gid_to_sat_gsl_if_idx,
        sat_neighbor_to_if,
        prev_fstate,
        enable_verbose_logs
):
    if enable_verbose_logs:
        print("  > Calculating Floyd-Warshall for graph including ground-station relays")
    ## Forwarding state
    fstate = {}
    num_node=num_satellites+num_ground_stations

    output_filename = output_dynamic_state_dir + "/fstate_" + str(time_since_epoch_ns) + ".txt"
    if enable_verbose_logs:
        print("  > Writing forwarding state to: " + output_filename)

    ## data structure convert
    adjacency_matrix=nx.to_numpy_array(sat_net_graph,nonedge=np.inf,dtype=np.float32).reshape(num_node**2)
    np.fill_diagonal(adjacency_matrix,0)
    dist_node_to_node_dev=floyd_in_parallel.floyd_gpu(adjacency_matrix)

    node_neighbor_index_list=np.cumsum(sum(adjacency_matrix!=np.inf))
    node_neighbor_array_list=np.array(np.where(adjacency_matrix!=np.inf))

    num_node_neighbor=node_neighbor_index_list[-1]
    num_int=num_node+num_node_neighbor+num_satellites**2+num_ground_stations+num_satellites+(num_node**2)*3
    num_float=num_node**2

    data_int=np.full(num_int,-1,dtype=np.int32)
    data_float=np.full(num_float,-1,dtype=np.float64)

    index_start=0
    for i in range(num_satellites):
        data_int[index_start]=node_neighbor_index_list[i]
        index_start=index_start+1

    for i in range(num_node_neighbor):
        data_int[index_start]=node_neighbor_array_list[1][i]
        index_start=index_start+1
    
    sat_to_sat_if=sat_neighbor_to_if.reshape(num_satellites**2)
    for i in range(num_satellites**2):
        data_int[index_start]=sat_to_sat_if[i]
        index_start=index_start+1
    
    for i in range(num_ground_stations):
        data_int[index_start]=gid_to_sat_gsl_if_idx[i]
        index_start=index_start+1

    for i in range(num_satellites):
        data_int[index_start]=num_isls_per_sat[i]
        index_start=index_start+1

    index_start=0
    for i in range(num_node**2):
        data_float[index_start]=adjacency_matrix[i]
        index_start=index_start+1

    ##copy data to device
    data_int_dev=cuda.to_device(data_int)
    data_float_dev=cuda.to_device(data_float)

    ## extract data
    index_start=0
    node_neighbor_index_dev=data_int_dev[index_start:index_start+num_node]
    index_start=index_start+num_node
    node_neighbor_id_dev=data_int_dev[index_start:index_start+num_node_neighbor]
    index_start=index_start+num_node_neighbor
    sat_to_sat_if_dev=data_int_dev[index_start:index_start+num_satellites**2].reshape(num_satellites,num_satellites)
    index_start=index_start+num_satellites**2
    gid_to_sat_gsl_if_dev=data_int_dev[index_start:index_start+num_ground_stations]
    index_start=index_start+num_ground_stations
    num_isls_per_sat_dev=data_int_dev[index_start:index_start+num_satellites]
    index_start=index_start+num_satellites
    next_hop_decision_matrix_dev=data_int_dev[index_start:index_start+(num_node**2)*3].reshape(num_node,num_node,3)

    index_start=0
    adjacency_matrix_dev=data_float_dev[index_start,index_start+num_node**2].reshape(num_node,num_node)
    ## invoke kernel function
    TPB=(32,32)
    row=int((num_node+31)/32)
    col=int((num_node+31)/32)
    BPG=(row,col)
    calculate_shortest_path_only_gsl_between_node[BPG,TPB](
        num_satellites,
        num_ground_stations,
        node_neighbor_index_dev,
        node_neighbor_id_dev,
        adjacency_matrix_dev,
        dist_node_to_node_dev,
        sat_to_sat_if_dev,
        gid_to_sat_gsl_if_dev,
        num_isls_per_sat_dev,
        next_hop_decision_matrix_dev
    )

    

    ## copy result to host
    next_hop_decision_matrix=next_hop_decision_matrix_dev.copy_to_host()

    with open(output_filename,"w+") as f_out:
        for src_node_id in range(num_node):
            for dst_node_id in range(num_node):
                next_hop_decision=(
                    next_hop_decision_matrix[src_node_id,dst_node_id][0],
                    next_hop_decision_matrix[src_node_id,dst_node_id][1],
                    next_hop_decision_matrix[src_node_id,dst_node_id][2]
                )

                if not prev_fstate or prev_fstate[(src_node_id,dst_node_id)]!=next_hop_decision:
                    f_out.write("%d,%d,%d,%d,%d\n" % (
                                src_node_id,
                                dst_node_id,
                                next_hop_decision[0],
                                next_hop_decision[1],
                                next_hop_decision[2]
                            ))
                
                fstate[(src_node_id,dst_node_id)]=next_hop_decision
    return fstate